Quantum entanglement of decohered two-mode squeezed states in absorbing and 

amplifying environment 
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We investigate the properties of quantum entanglement of two-mode squeezed states interacting 
with linear baths with general gain and loss parameters. By explicitly solving for p from the master 
equation, we determine analytical expressions of eigenvalues and eigenvectors of p TA (the partial 
transposition of density matrix p). In Fock space, p TA is shown to maintain a block diagonal 
structure as the system evolves. In addition, we discover that the decoherence induced by the baths 
would break the degeneracy of p TA , and leads to a novel set of eigenvectors for the construction of 
entanglement witness operators. Such eigenvectors are shown to be time-independent, which is a 
signature of robust entanglement of two-mode squeezed states in the presence of noise. 

PACS numbers: 03.67.Mn, 03.65.Yz, 42.50.Dv, 42.50.Lc 



I. INTRODUCTION 

Optical two-mode squeezed vacuum (TMSV) has been 
a major source of continuous-variable entanglement for 
quantum communication In recent years, intriguing 
applications such as quantum teleportation 0, H, 0,TE Q 
and quantum dense coding [3, HJ have been demon- 
strated experimentally with TMSV. Theoretically, it is 
also known that TMSV maximizes the EPR correlation 
when a fixed amount of entanglement is given [?|. In 
order to exploit fully the non-classical properties of such 
entangled light fields, it is important to understand deco- 
herence effects as they propagate through noisy environ- 
ments US El GJ, [H El EE El EE GJ, El • This belongs 
to a more subtle topic involving the characterization and 
quantification of mixed state entanglement in general. 

For bipartite systems, Peres and Horodecki have de- 
veloped a powerful criterion of entanglement, which is 
known as the PPT (positive partial transposition) crite- 
rion 0, [2l], [2^]. If the partial transposition of a den- 
sity matrix (denoted by p Ta ) has one or more negative 
eigenvalues, then the state is an entangled state. Phys- 
ically, the partial transposition for separable states can 
be considered as a time-reversal operation, and one can 
construct a variety of uncertainty relations serving as in- 
dicators of entanglement [H [2J] . For two- mode Gaus- 
sian states such as TMSV, PPT provides a necessary and 
sufficient condition of separability [H [26[ . 

The dynamics of disentanglement of TMSV in various 
noisy situations has been addressed byseveral authors 
recently [H El El El, El, El El EE El, M- Th e 
fact that an amplitude damped TMSV remains Gaus- 
sian enables an elegant description of entanglement based 
on the properties of covariance matrix associated with 
the density operators [27| . In particular, from the time- 
dependent solution of Wigner function fl5|. 16l. E3] or 
the corresponding characteristic function 
quantify the degradation of entanglement b 
the negativity [28j and relative entropy [2J 
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one can 



calculating 
. It is now 

known that for an initial TMSV at a non-zero thermal 



bath, quantum entanglement vanishes completely in a fi- 
nite time [27j | . 

However, we notice that there are much less investiga- 
tions addressing the structure of p Ta directly, and yet p TA 
is what the PPT criterion originally based upon. Since 
P Ta could manifest differently in various basis, the study 
of p TA in Fock space, for example, could reveal entangle- 
ment properties not easily found by the Wigner function 
method [30( . An example we notice is the construction of 
entanglement witness operators via the projectors formed 
by the eigenvectors of p TA with negative eigenvalues [3l[ . 
Such entanglement witness operators, which correspond 
a variety of observables for the detection of entanglement, 
are determined by p TA . 

The main purpose of this paper to indicate some key 
features of decohered entanglement as revealed by eigen- 
values and eigenvectors of p Ta . Our analysis will con- 
centrate on the structures of eigenvectors in Fock space, 
which is also where interesting non-local correlations of 
continuous- variable systems can be observed 32]. For 
a TMSV under the influence of amplitude damping (or 
gaining in an amplifier), we solve for the time evolution of 
p and determine the exact eigenvectors and eigenvalues 
of p TA analytically. These eigenvectors are shown to have 
a strong correlation in photon numbers, and hence p TA is 
a block diagonal matrix in Fock space. Therefore witness 
operators associated with each block involve only a finite 
number of Fock vectors, which implies that the detection 
of entanglement can only require a small portion of the 
Hilbert space. This is in contrast to entanglement detec- 
tion based on uncertainty relations in which the entire 
Hilbert space is usually involved [H H, El HI • In this 
sense the eigenvectors of p TA access the entanglement sig- 
natures 'locally', which is a complement to 'global' char- 
acterization (of Gaussian states) using covariance matri- 
ces. As we shall see below, as long as the initial state is 
a TMSV, the corresponding eigenvectors do not change 
with time, indicating that the entanglement carried by 
TMSV is robust against amplitude damping. 
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FIG. 1: The sub-matrix structure of p T + with initial TMSV, 
in the Fock basis. 



II. MASTER EQUATION AND SOLUTION 

To begin with, we consider the time evolution of 
an initial TMSV, each coupled with a separate phase- 
insensitive linear bath. In terms of the annihilation op- 
erators a and b of the two modes, the master equation 
governing the dynamical process is [33| : 

p = G(2a) pa — <W p — paa) + 26' pb — bw p — pbb^) 



L(2apa) — a) ap — pa) a 



2bptf -tfbp- ptfb), 



(1) 



where G and L are the gain and loss parameters respec- 
tively, both having a dimension of time -1 . Depending 
on the values of G and L, the master equation describes 
amplifying or damping effects due to the coupling with 
the baths. For dissipation in thermal baths, each of 
temperature T, we have the parameters G = ^n t h and 



(n t h + 1), where n th 



— is the average 



2 v-tn i ,„ m — exp ( ha! / kT y 

number of photons in each of the modes (with frequency 
lj) at thermal equilibrium, and 7/2 is the decay rate of 
the mode amplitudes. In this paper we focus on the ini- 
tial TMSV with the squeezing parameter r > 0: 



|V(0)) = exp[r(aW - ah)] |00) = \fl - A 2 ^ A™ \nn) ,(2) 

n=0 

where A = tanhr and |00) is the two-mode vacuum state. 

A. Block structures of p TA in Fock space 

To investigate the entanglement properties of the den- 
sity matrix p, we study its partial transposition p TA . The 
p TA is of infinite dimension, however, by examining the 
master equation in the Fock basis, block structures of 
p TA can be identified. Let us denote the matrix elements 
of p TA by 



P Ta 



(nm\p TA \pg) — (pm\p\nq} , 



which is governed by the following differential equation: 

+2^fmqp TA nm _ lpq _ 1 
-(n + m + p + q + A)p TA nmpq ] 
+L[2^(n + l)(p + l)p TA n+hm , p+hq 
+2^(m + l){q + l)p TA 
~(n + m + p + q)p TA nmpq \. (4) 



and 



P 



n,m,p,q 



(t = 0) = 5 pm S nq (l - A 2 )A 



(5) 



corresponds to the initial condition (2). 

It can be seen from Eq. |4} that each el- 
ement P TA nmpq(t) ^ s coupled with elements 
p TA n +l,m+k,p+l,q+k(°) onl y> for integers I and k. 
Therefore P Ta n+l , m+k , m +Ln+k (*) arc thc onl Y non-zero 
elements at any time t > because of the initial con- 
dition. By noting that the sum of the first two indices 
equal to that of the last two, we can group all non-zero 
elements P Ta n+ i, m+k , m+ i, n+k (t) into sub-matrices M s 
according to the sum index S = n + / + m + k, i.e., 
P TA n +l,m+k,m+l,n+k (<) is contained in M n+ i +m +k- We 
can therefore express p TA in a direct sum of Ms a s 
follows: 



(6) 



5=0 



where the sub-matrix Ms has a dimension of 5+1, since 
elements in Ms have its first two indices as {0, S}, {1, 5— 
1}, ... , {5,0}. Fig. [T]shows the sub-matrix structure of 
p TA . Note that characteristic sum S is equal to the total 
number of photons that the two modes contain. From 
Eq. (|U) we observe that probabilistic flow occurs between 
elements in neighboring sub-matrices, with emission or 
absorption of one photon in one of the modes at one 
time. 

The time evolution of a typical sub-matrix of p TA is 
illustrated schematically in Fig. [2l which will be dis- 
cussed in detail in the later part of the paper. Ini- 
tially, only opposite-diagonal elements are present, hav- 
ing the magnitude as A s . As time increases, element 
flows from neighboring sub-matrices, and disentangle- 
ment of the sub-matrices occurs at a critical time t = 
t c (Section IIII A[) . In the case of thermal bath, p TA 
evolves into a diagonal form in the long time limit, 
settling as the thermal equilibrium state p TA nmpq = 

A A 1 I nth \n+m 

OnpOmq (n t(l + l) 2 v n th + l > 



(3) 



B. Analytic solution of p in position space 

To analyze the properties of p TA , it is more conve- 
nient to first determine p in position space and then 



3 





l = Q 



CO 



FIG. 2: (Color online) Schematic diagram showing the evo- 
lution of the distribution of elements in a sub-matrix of p TA , 
assuming the baths are thermal baths. The t c is the critical 
time for disentanglement. 



make the transformation to Fock space. The position 
space method of finding p was previously employed in 
Ref. [ll|, H3 m studying entanglement in various oscil- 
lator systems. In this subsection, we present an explicit 
solution of master equation (1) with an initial TMSV. 
We remark that our method is different from that given 
in [ll[, as the latter involves a Fourier transform of the 
density matrix, i.e., the momentum space. Here we solve 
the density matrix entirely in position space (Appendix 
A). This turns out to be more convenient for the real 
symmetric Gaussian states considered here, since fewer 
differential equations are involved. In addition, the re- 
sultant solution is more transparent for further analysis 
of eigenvectors in the next section. 

Let us denote the 'position' operators as x = -^(a+a^) 

and y — + and define 

p(xi,yi;x 2 ,y 2 ;t) = (xi,yi\p(t) \x 2 ,y2) , (7) 
then the master equation (1) becomes, 



P = -^\_L{x\ + x\- 2xix 2 +yi + y£- ^ViV2 - 4 



• rr - d z — cT -2(9 <9 

X2 VI V2 AU Xl U X 2 



28 8 

6U yi Uy 2 



-2xid X2 - 2x 2 d Xl - 2yid y2 - 2y 2 d yi ) 
+G(xf + x 2 2 - 2x\x 2 + y\ + y 2 2 - 2y x y 2 + 4 
-fi 2 _ f) 2 - f) 2 - f) 2 _ 2<9 B - 2(9 (9 

u xi u x 2 u yi u y 2 Au yi u V2 

+2 Xl d X2 + 2x 2 Xl + 2 Vl d V2 + 2y 2 d yi )]p. (8) 

For an initial state (2), p{x\,yi, x 2l y 2 ; i) takes a Gaussian 
form at any time t, 

p{xi,yi;x 2 ,y 2 ;t) = E(i) exp[-A(t)(xl + x% + y\ + yj) 
+B(t)(x 1 yi + x 2 y 2 ) 
+C(t)(xix 2 + yiy 2 ) 
+D{t){x lV2 + x 2Vl )l (9) 

where A(t), B(t), C(t) and D{t) are real time-dependent 
coefficients, and the normalization factor is: 



S(i) = ^[2A{t)-C{t)] 2 -[B(t)+D{tW 



(10) 



By substituting Eq. ^ into the master equation, the 
coefficients are found to obey a set of coupled equations 



that can be solved analytically (Appendix [A} . For the 
TMSV considered here, we have, 

1 C A- J (x 2 ) 

Ait) = -\2A r] + — — (n-l) + 5-]. 

U 4 [ G-L [1 > 2((x 2 ) 2 - (xy) 2 ) h 

r 2((x 2 ) 2 - (xy} 2 t ) h 

1 C 4- T (x 2 ) 

C(t) = -[2A oV + ^-—(v-l) — rl. 

2 G - L 2{{x 2 ) 2 - {xy) 2 t ) b 



2 [ 2((x 2 ) 2 - (xy) 2 ) 1 



(11) 



it \^y/t' 

Here A$ = | cosh2r, Bq = sinh2r and 77(f) = exp[2(G — 
L)t] are defined, and the expectation values are given by, 

(x 2 ) t = A t 1 + n G ^ +L r , 



* ' 2{G-L) 

(xy)t = -y 7 /- 

III. PROPERTIES OF P Ta 



(12) 



According to PPT criterion, the appearance of nega- 
tive eigenvalues of p TA is a signature of entanglement. In 
this section, we solve the eigenvectors and eigenvalues of 
P Ta as the system evolves. Then we discuss how deco- 
herence affects the entanglement properties of p Ta . The 
eigenvalues and eigenvectors of p TA are defined by: 

P Ta (xi 1 y 1 ;x 2l y 2 ;t) ip rhm (x 2 , y 2 ; t)dx 2 dy 2 

= Cn,m(t)Vn,m(xi,yi;t). (13) 

Our main technique of solving the eigen-problem is the 
use of Mehler formula which expands a double Gaussian 
function into a series of orthogonal functions. After some 
calculations (see Appendix |B|) . we obtain the expression 
of eigenvalues, 



€n,m(t) = E(t)7T- 



and the corresponding eigenvectors, 



r (i4) 



1 



V2"+ m n!m!7r 



Xl 



V2 



Vl )H m ( Xl 



yi 



V2 



1 



xexp[-^(xl+yf)] 



(15) 



where H n are the Hermite polynomials. In writing 
Eq. HH), we have defined 



ai(t) = 




- Bit) + 


C(t) 


+ D(t)}, 


Pl(t) = 




+ Bit) - 


C(t) 


+ D(t)}, 


a 2 (t) = 




+ B(t) 4 


C(t) 


-D{t)l 


Mt) = 




-B(t)- 


C(t) 


-D(t)} (16) 
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and the solution of A, B, C and D are given by Eq. (fTTj) . 
Note that in writing Eq. (fT5|) from (B8), we have used 
the fact that for TMSV, aifli = CC2P2 = Jq for all time 
t > 0. 

We now transform the eigenvectors from the position 
space to the Fock space. Note that p TA is a basis de- 
pendent operation. The eigenvectors of p TA defined in 
two different basis sets do not transform directly. An 
exception is the case when the two sets of basis vectors 
transform by a real unitary matrix [35j ] , which is the case 
here. This allows us to write down the eigenvector |<^n,m) 
in Fock space from Eq. (TTF 



\Yr, 



Vn! V2 



V V2 



-moo), a?) 



which can be connected to the sub-matrices of p TA (Fig. 
[T]) via the photon number sum S = m + n to label the 
eigenket, so that 



1 



\<Pn,S- 



y/¥^/n\{S-n)\ 

s 

x £ VMS-jV- \j, S-j). (18) 



Here we have used the abbreviation 



T S,n,j 



min(j,r 

E 

k=0 



(_ 1) n-fc Cr S-n cs?j 



(19) 



with = 



si 



In this way \ip n ,S-n) an d £ n> , 



rl(S-r)! 

(n = 0, 1, S) are the n th eigenvector and eigenvalue of 
the block with characteristic sum S. 



A. Evolution of negative eigenvalues 

By inspecting Eq. (fl4")) . we find that all the negative 
eigenvalues in each block share the same value at t = 
(Fig. [3]). The same is true also for positive eigenval- 
ues. However, such a strong degeneracy is broken by 
coupling with the baths. This is illustrated in Fig. [3] 
where the time-dependence of individual eigenvalues in 
various block indices S is shown. Except at t = and at 
the critical time t c , we see that the negative eigenvalues 
possessing different values. 

It is important to observe that the eigenvalues are neg- 
ative for odd n, when y/ai < y/Jh, or in other words, 
B(t) > C(t). All eigenvalues turn zero at the same time, 
when we have B(t) = C(t), except for the only eigenvalue 
with n = in each block. Such a critical time t c is given 
by, 



tc = 



1 



G + LX 



2{L-G) a \G{l + X) 



■log 



(20) 



which is always positive finite as long as G ^ 0. For 
the case G — > 0, we have t c — > 00. We remark that the 
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FIG. 3: (Color online) The eigenvalues of sub-matrices of p TA 
of an amplifier system with G = 1.57, L = O.57 and initial 
squeezing factor r = tanh" 1 0.2, with characteristic sum (a) 
S = 1, (b) S = 4, (c) 3 = 7, and (d) S = 10. Red line 
indicates the n = eigenvalue which does not turn zero at 
t c , while blue line indicates the most negative eigenvalue with 

71= 1. 



disentanglement time t c was previous obtained in Ref. 
[26| for thermal baths, here we obtained a general ex- 
pression (|20p that applies to linear amplifiers as well. In 
particular, in the case when gain and loss parameters are 
equal, i.e., G = L, the critical time can be reduced to 

f — 1 A > f) 
Lc ~ 2G l+A ^ u - 

We point out that at the time of disentanglement t = 
t c , there is only one non-zero eigenvalue (with the index 
n = 0) in each sub-matrix (Fig. [3]). Therefore p TA at the 
critical time is highly degenerate, and the corresponding 
symmetry property of p TA (t — t c ) is indicated in the 
relation: 



n TA 

P j,S-j,j,S-j 



n TA 

P j,S-j,S-j,j 



~P S-j,j,j,S-j - P S-j,j,S-j,j- 



(21) 



This results in the symmetric distribution of elements as 
shown in Fig. fj] schematically. 

As a further remark, it is interesting that negative 
eigenvalues may not necessarily be monotones over time. 
This can be seen by differentiating Eq. (fT4"|) and looking 
at the initial rate: 



t=o 



(—1)™ tanh s r csch rsech 3 r 

x{(G- L)(S -2n) 
+ (G + L)(S-2n) cosh2r 
-[LS + G{2 + S)]smh2r}, (22) 



which can result in a negative value for certain param- 
eters, i.e., some eigenvalues of odd n can become more 
negative over time (Fig. [5Ji). An exceptional case is when 
G = 0, where we find that for odd n, the derivative at 
t = must be positive by inspecting Eq. 
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FIG. 4: (Color online) The negativity Af of p TA with different 
G, fixing the parameters L — 7 and initial squeezing factor 
r = tanh" 1 0.2. 



B. Negativity and sub-negativity 

Negativity Af serves as a computable measure of en- 
tanglement denned by the trace norm of p TA minus 1 
divided by 2 (28|. For separable states, p TA is still a den- 
sity matrix with trace 1 and hence Af — 0. However, for 
non-separable states with negative p TA , we have Af > 0. 
Specifically, Af equals the sum of the absolute value of 
negative eigenvalues of p TA [H, [36| . For the system con- 
sidered in this paper, Af reads, 



AT 



wE(t) 



1 



1 



(23) 



Alternatively, Af can be derived from the symplectic spec- 
trum of the covariance matrix associated with the den- 
sity operator [27| ■ In Fig. [4] we show the time evolu- 
tion of negativity Af for initial TMSV with different G 
parameters. An example for the G — case is the zero- 
temperature bath dissipation scenario. Fixing L, we ob- 
serve that a larger gain G leads to a shorter t c . 

We can also calculate the negativity in a sub-matrix 
A4s, which measures the contribution of entanglement 
from the corresponding block that builds up p TA . Specif- 
ically, the sub-matrix negativity Afs is defined the same 
way as negativity Af but restricted to the sub-matrix of 
p TA with the characteristic sum S. To our knowledge, 
such a sub-matrix negativity, which requires the calcu- 
lations of individual eigenvalues, has not been discussed 
before. Explicitly, Afs takes the form: 



\2n+l 



S-2n-l 



(v / aT+ VPl) 



2n+2 (V^ + VJh) s ~ 2n 



(24) 

where P is defined as the integral part of In Fig. 

[5] the behavior of Afs of some blocks is illustrated. It 
is surprising that that for higher sub-matrices the corre- 
sponding negativity Afs may increase over time. Physi- 
cally, the increase of Afs is due to the probability flow in 




FIG. 5: The negativity of sub-matrices Afs of p TA , for an am- 
plifier system with G = I.57, L — O.57 and initial squeezing 
factor r — tanh -1 0.2, with characteristic sum (a) 5=1, (b) 
S = 4, (c) S = 7, and (d) S = 10. 



Fock space arising from damping or amplifying mecha- 
nisms. Since each sub-matrix of p Tjl does not necessarily 
conserve probability (i.e., the trace of A4s is not a con- 
stant), it is possible that some blocks could have their 
negativity increasing with time. 

However, the increase of Afs does not violate the fact 
that the overall negativity Af of p TA is an entanglement 
monotone that does not increase under LOCC (the mas- 
ter equation corresponds to local operations). Eq. (fl4|) 
reveals that eigenvalues of higher blocks are of smaller 
order of magnitude. As we see in Fig. O although nega- 
tivity of individual higher sub-matrices may increase over 
time, their contribution for negativity is smaller by sev- 
eral orders than the sub-matrices with lower S, and there- 
fore the overall negativity is still monotonic decreasing. 



C. Robust structure of entanglement witness 

An entanglement witness operator W is designed for 
the detection of entanglement such that Tr(pW) < 
for some non-separable states p, but Tr(p sep W) > 
for all separable states p sep [2l|. For each eigenvec- 
tor \4>) of p TA with a negative eigenvalue, one can con- 



struct an W by W = \<f>)(cf)\ TA [3l| meeting the criteria 
above. In our system, we can construct a family of W 
from eigenvectors \(p n ,s-n) with odd n accordingly, i.e., 
Ws> = \<Pn,S-n) {<fn,s-n\ TA - From Eq. (Tn§, the explicit 
form of Ws.n reads, 



s s 



1 



2 s nl(S-n)l 



Fs,n,jFs,n,l 



j=0 1=0 

X y/j\(S-j)\l\(S-l)\} \l, S - j) (j, S-l\, 

(25) 
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FIG. 6: Example of entanglement witness W2,i- 



which shows that Ws,n operates in only a finite dimen- 
sion in the Fock space. Note that Tr(pWs, n ) — £n,s-n, 
and is the most negative eigenvalue in each sub- 

matrix, therefore Ws.i provides the most significant en- 
tanglement detection among all witnesses constructed 
from vectors lying within the sub-matrix Ms- Some ex- 
amples for entanglement witnesses with n = 1 are shown 
in Fig. [S]and Fig. El 

We observe that W also has a block diagonal structure 
in which non-zero elements are: (j, S — l \ YVs.n \l,S — j). 
This allows us to divide W into sub-blocks Ak , with each 
sub-block characterized by a difference Kg [— S, S]: 



K=-S 



K- 



(26) 



where the element (j, S — l\ Ws, n \l, S — j) lies in the sub- 
block with K = S — j — I. The explicit form of Ape is 



A 



K 



min(S,S-K) 

E { 

;=max(0,-.R") 



1 r r 

z^n'^S — to)! 



x ^{S -K- l)l(K + l)W(S - /)!} 
\l,K + l){S-K-l,S-l\, (27) 

having a dimension of S — \K\ + 1. We remark that the 
simplest 2x2 entanglement witness Wi,i was constructed 
in [37J using a different approach. Here our general Ws, n 
applies to all S and n. 

Finally, let us emphasize the robustness feature of 
TMSV against decoherence. We have seen from Eq. (fTS)) 
explicitly that the eigenvectors \<p n ,s~n) remain un- 
changed with time. The time-independent \(p n .s-n) sug- 
gests that TMSV is robust against noise, in the sense that 
structure of entanglement witness Ws, n is preserved. The 
degradation of entanglement would only affect the eigen- 
values. We stress that such a time-independent property 
of eigenvectors is specific to initial TMSV, and does not 
hold for arbitrary initial states in general. In Appendix 
IB1 we derive the eigenvectors of p TA evolving from an ini- 
tial two-mode symmetric Gaussian states with arbitrary 



FIG. 7: Part of the entanglement witness Ws t i, showing the 
blocks with characteristic difference from -3 to 0. W3.1 is 
symmetric about the top-right to bottom-left diagonal. 



real coefficients Aq, Bo, Cq and Do- We find that the 
time dependence of the eigenvectors arises solely from 
the evolution of the factors a\(3\ and a2^2- Because of 
the system-bath interactions, these factors are generally 
dependent on time. However, the state evolving from 
TMSV is an exception in which one can show that the 
corresponding a\f3\ and a2/?2 equal the constant 1/16 at 
all times. 



IV. CONCLUSION 

To summarize, by deriving an exact analytic solution 
for p TA and examining its eigenvectors, we discover sev- 
eral important features about the loss of entanglement 
of a TMSV suffered from decoherence. Both amplitude 
damping and amplification effects have been included in 
our analysis. Throughout the decoherence process, the 
block diagonal structure of p TA is shown to be main- 
tained in Fock basis. As each block spans only a finite 
portion of the Fock space, the existence of negative eigen- 
values in a block implies that entanglement information 
can 'survive' in the corresponding photon-number sub- 
space. If quantum entanglement of the system is to be 
destroyed completely, then all the blocks have to be made 
positive. In other words, by simply mixing the system 
with another state involving finite Fock space would not 
destroy the entanglement. For the decoherence process 
considered in this paper, all the blocks turn positive at 
a critical time t c , which agree with the previous analy- 
sis based on the covariance matrix. At t < t c , we de- 
rived an explicit expression of the negativity as well as 
the negativity of sub-matrices in order to characterize 
the time-dependence of entanglement. Interestingly, the 
negativity of some sub-matrices could increase with time 
when G is non-zero, although the effect is weak accord- 
ing to our calculations. The signature of entanglement in 
photon number subspace can be detected by the witness 
operators Ws. n , and we have constructed Ws,n explicitly 



7 



in this paper. These witness operators are also block di- 
agonal, and most remarkably, they are time-independent 
even though the system is under the influence of noise. 
We interpret such a property as a kind of robust entan- 
glement structure inherited in TMSV. 
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APPENDIX A: TIME EVOLUTION OF 
GENERAL REAL SYMMETRIC TWO-MODE 
GAUSSIAN DENSITY OPERATOR 

In this Appendix, we consider the time evolution of 
real symmetric two-mode Gaussian states, each coupling 
linearly with a separate bath with general gain and loss 
parameters. The density matrix is represented in position 
space, obeying the master equation as in Eq. ||5J). For the 
solution of the two-mode real symmetric Gaussian state 
presented in Eq. ©, a direct substitution leads to the 
following coupled differential equations: 



.4 

B 
C 
D 



1 



(G - L)C + ^(G + L)[1- (2A - Cf 



(B + D) 2 }, 

-2(G - L)D - 2(G + L)(2A - C)(B + D). 

4(G - L)A + (G + L)[1 + (2A - Cf + (B + D)% 

-2(G-L)B-2(G + L)(2A-C)(B + D). (Al) 



Without loss of generality, we consider the G ^ L case. 
From Eqs. (Al), B - D and 2 A + C have the simple 
solution, 

B(t)-D(t) = (Bo-Do)v(t), 

2A(t) + C(t) = (2A + CoMt) + ^±4[ v (t) - l], 



G 



(A2) 



where rj(t) = exp[2(G — L)t] and the zero subscripts de- 
note the values at t = 0. The other two combinations, 
B + D and 2 A — C, can be found by noting that they are 
related to the second moments (a; 2 ) and (xy) by: 



2A(t)-C(t) = 
B(t)+D(t) = 



2{{x 2 ^ 



• t ( x y)t) 
(xy) t 



2\(x 2 



(A3) 



with the subscript t denoting the value at time t. Under 
the condition: 



(« 2 ) t = (*) t = (*% = 



(A4) 



which applies to TMSV, and from the Heisenberg equa- 
tions of motions of a and b, the time-dependence of the 



second moments are given by: 

(xy) t = {xy) Q r,{t). (A5) 
Here the initial second moments are given by: 

2A - Co 



<* 2 >o = 

(xy) = 



2[(2A -Co) 2 -(B + D yy 

Bq + Dq 

2[(2A - C ) 2 - (B + D ) 2 }' 



(A6) 



Thus the time evolution of the coefficients of the two- 
mode Gaussian state solution are as follows: 

A(t) = I[(JL4o + a,)»7+§3^fo-l) 



<* 2 > t 



2((x 2 ) 2 - (xy) 2 )* 



2 2((x 2 ) 2 - (xy) 2 ) 1 

C(t) = l[(2Ao + C )r,+ ^±^(r,-l) 



2((x 2 ) 2 -(xy) 2 y , 



D(t) = h-(Bo-D )v+ ^ t-], (A7) 

2 2((x 2 ) 2 - (xy) 2 y 



where the expectation values are given in Eq. (|A5[) . 

In the case of TMSV, the solution is reduced to Eq. 
(12) as the initial coefficients satisfy: 



Go = Do = 
AA 2 - B 2 = 1. 



(A8) 
(A9) 



In particular, by noting that Ao = (x 2 ) Q = \ + ( a ^ a ) = 
h + r^ra , we have A — \ cosh 2r. 



APPENDIX B: DERIVATION OF 
EIGENVECTORS AND EIGENVALUES OF P Ta 

In position space, p TA (xi,yi;x 2 ,y 2 ;t) = 
p(%2, yi',Xx,y2',t). The eigenvectors tp n m and eigenvalues 
£,nm of p TA are defined by Eq. (fl3|) . Applying the trans- 
formation Xj — (uj + Vj)/\f2 and y.j = (—uj+ Vj)/y2 
(j = 1,2), p TA becomes a neat product of two double 
Gaussians as follows: 

p TA (u 1 ,vi;u 2 ,v 2 ;t) 

= E(t) exp[-ai(i)(?ii - u 2 ) 2 - Pi{t)( Ul + u 2 ) 2 } 
x exp[-a 2 (t)(ui - v 2 f - (3 2 (t){ Vl + v 2 )% (Bl) 

where ai = \{2A- B + C + D), f3 x = {{2A + B-C + D), 
a 2 = \{2A + B + C - D) and j3 2 = \{2A — B — C — D). 
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We apply Mehler's Formula twice, one for the u\, u 2 
double Gaussian, and one for the V\ , v 2 double Gaussian 
function. This would lead to the Schmidt decomposition 
on p 



Ta- 



p TA (u 1 ,v 1 ;u 2 ,v 2 ;t) = S(i)-(- 



-r 



4 ai/3ia2/?2' 

oo oo 
X X! ^nfn{u\)fn{u 2 ) ^ fm(vi) f m (v 2 ) , (B2) 
n— m— 

where the Schmidt modes f n (u) and f m (u) are: 
1 



/n(") 



fm(u) 



(- )*H„[2(ai/3i)4 U ] 



\ / 2 n - 1 n\ 7T 
exp(— 2y ai/3iw ) 



' :(^)'i/ m [2(^)^] 



V2™ rT m! 7T 



exp(-2 v /a 2 /?2U 2 ), 
and the coefficients A„(t) and A m (t) 



(B3) 



A„(t) = 2( ai /3i)4 



A m (t) =2(a 2 /3 2 )3 



(V^T+^T) n+1 ' 



(V^ + V&) m+1 ' 



(B4) 



where H n (u) are the Hermite polynomials. Rearranging 
terms, Eq. (|B2|) gives: 



p TA (xi,yi;x 2 ,y 2 ;t) = £ w , m ft) y n,m (ad , 2/1 ; *) 

(x 2 ,y 2 ;t), (B5) 



where the eigenvectors <^ n ,m are 



ip n ,m{zi,yi;t) 
1 



V2"- 1 2 m - 1 n!m! 



( =2 ) 4 



x exp[- y/aifii (xi - yi) 2 ] exp[-v/a 2 /3 2 (a;i + yi) 2 ], 

(B6) 



and the eigenvalues £„ jm of p TA are, 



Cn,m(t) = S(t)7T 



m+1 



(B7) 



These expressions of eigenvectors and eigenvalues are 
for any time-dependent coefficients A(i), B(t), C(t) and 
D(t), i.e., applicable to states evolving from arbitrary 
initial values Aq, Bq, Cq and Dq. The special case with 
initial TMSV is given in Eq. flU]) and Eq. (fig) . 
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